1. Data Preparations

library(zoo)
data <- read.table("BTC-USD-2013-2017.csv", sep=",", header=TRUE); attach(data) # original data
dates <- seq(as.Date('2013-01-01'), by = 'days', length = 1826) # transform to date formats
data$Date <- dates # append to original data
data <- read.zoo(data) # zoo formats
plot.zoo(data[,1], main = "BTC/USD", xlab = "Year", ylab = "Dollars")

plot.zoo(data, main = "all variables", xlab = "Year") 

2. Bitcoin Data

2.1 Visualize Data

library(astsa);library(PerformanceAnalytics); library(evir);library(zoo)
bitcoin.logreturn <- diff(log(data[,1]))*100
bitcoin.logreturn.squared <- bitcoin.logreturn**2
par(mfrow=c(1,2))
plot.zoo(data[,1], main = "BTC/USD", xlab = "Year", ylab = "Dollars")
plot.zoo(bitcoin.logreturn, main = "BTC/USD Log Return", xlab = "Year")

par(mfrow=c(1,1))
acf(coredata(bitcoin.logreturn),lag.max = 200)

acf(coredata(bitcoin.logreturn**2), lag.max = 100);pacf(bitcoin.logreturn)

bitcoin.logreturn.squared <- bitcoin.logreturn**2

2.2 Check Normality

par(mfrow=c(1,2))
qqnorm(bitcoin.logreturn);qqline(bitcoin.logreturn)
hist(bitcoin.logreturn,main = "Histogram", xlab = "Return")

2.3 Fit an ARMA Model Using auto.arima

library(forecast)
auto.arima(bitcoin.logreturn, seasonal = FALSE, ic = "bic")
## Series: bitcoin.logreturn 
## ARIMA(3,0,2) with zero mean 
## 
## Coefficients:
##          ar1      ar2     ar3      ma1     ma2
##       0.8988  -0.8073  0.2676  -0.8834  0.5792
## s.e.  0.0544   0.0419  0.0261   0.0544  0.0484
## 
## sigma^2 estimated as 40.57:  log likelihood=-5966.31
## AIC=11944.63   AICc=11944.67   BIC=11977.68

ARMA model: \[X_t = \]

2.4 Check Residuals of ARMA Model

library(astsa)
btc.logreturn.return <- sarima(bitcoin.logreturn, 3, 0, 2)
## initial  value 1.913161 
## iter   2 value 1.891052
## iter   3 value 1.880737
## iter   4 value 1.879690
## iter   5 value 1.877509
## iter   6 value 1.864927
## iter   7 value 1.863816
## iter   8 value 1.862380
## iter   9 value 1.862181
## iter  10 value 1.861953
## iter  11 value 1.861903
## iter  12 value 1.861815
## iter  13 value 1.861609
## iter  14 value 1.859861
## iter  15 value 1.858589
## iter  16 value 1.857562
## iter  17 value 1.857466
## iter  18 value 1.857038
## iter  19 value 1.853649
## iter  20 value 1.851540
## iter  21 value 1.850974
## iter  22 value 1.850466
## iter  23 value 1.849721
## iter  24 value 1.849507
## iter  25 value 1.849467
## iter  26 value 1.849459
## iter  27 value 1.849455
## iter  28 value 1.849454
## iter  29 value 1.849451
## iter  30 value 1.849450
## iter  31 value 1.849449
## iter  32 value 1.849449
## iter  33 value 1.849449
## iter  34 value 1.849449
## iter  35 value 1.849448
## iter  36 value 1.849448
## iter  36 value 1.849448
## iter  36 value 1.849448
## final  value 1.849448 
## converged
## initial  value 1.848747 
## iter   2 value 1.848747
## iter   3 value 1.848747
## iter   3 value 1.848747
## iter   3 value 1.848747
## final  value 1.848747 
## converged

btc.logreturn.return
## $fit
## 
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, 
##     Q), period = S), xreg = xmean, include.mean = FALSE, optim.control = list(trace = trc, 
##     REPORT = 1, reltol = tol))
## 
## Coefficients:
##          ar1      ar2     ar3      ma1     ma2   xmean
##       0.8885  -0.8059  0.2622  -0.8765  0.5784  0.3792
## s.e.  0.0555   0.0418  0.0264   0.0554  0.0486  0.1592
## 
## sigma^2 estimated as 40.34:  log likelihood = -5963.53,  aic = 11941.05
## 
## $degrees_of_freedom
## [1] 1819
## 
## $ttable
##       Estimate     SE  t.value p.value
## ar1     0.8885 0.0555  16.0191  0.0000
## ar2    -0.8059 0.0418 -19.2984  0.0000
## ar3     0.2622 0.0264   9.9299  0.0000
## ma1    -0.8765 0.0554 -15.8095  0.0000
## ma2     0.5784 0.0486  11.9080  0.0000
## xmean   0.3792 0.1592   2.3813  0.0174
## 
## $AIC
## [1] 4.70384
## 
## $AICc
## [1] 4.70497
## 
## $BIC
## [1] 3.721953

2.5 check residuals

2.5.1 Visualize Residuals and Squared Residuals

res.btc.logreturn.return <- resid(btc.logreturn.return$fit)
par(mfrow=c(2,2))
qqnorm(res.btc.logreturn.return);qqline(res.btc.logreturn.return);
hist(res.btc.logreturn.return); acf(res.btc.logreturn.return)
mtext("Residuals", side = 3, line =2, outer = TRUE)
res.square.btc.logreturn.return <- res.btc.logreturn.return**2
par(mfrow=c(1,2))

qqnorm(res.square.btc.logreturn.return);qqline(res.square.btc.logreturn.return);
acf(res.square.btc.logreturn.return,lag.max = 100)

#mtext("Squared Residuals ", side = 3, line =2, outer = TRUE)

2.5.2 ARCH Engle’s Test

library(aTSA)
arch.test(arima(bitcoin.logreturn, order = c(3,0,2)))
## ARCH heteroscedasticity test for residuals 
## alternative: heteroscedastic 
## 
## Portmanteau-Q test: 
##      order  PQ p.value
## [1,]     4 211       0
## [2,]     8 446       0
## [3,]    12 450       0
## [4,]    16 451       0
## [5,]    20 452       0
## [6,]    24 452       0
## Lagrange-Multiplier test: 
##      order   LM p.value
## [1,]     4 7952       0
## [2,]     8  999       0
## [3,]    12  570       0
## [4,]    16  399       0
## [5,]    20  309       0
## [6,]    24  253       0

Therefore, squared residuals are heteroscedastic.

Fit Garch Model

library(fGarch)
source("garchAuto.R")
#fit = garchAuto(bitcoin.logreturn, cores=8, trace=TRUE) ## can only work in MAC

btc.garch1 = garchFit(~arma(3,2)+garch(1,1), cond.dist = "std", data = bitcoin.logreturn, trace = FALSE, include.mean = FALSE)
btc.garch2 = garchFit(~arma(3,2)+garch(1,1), cond.dist = "sstd", data = bitcoin.logreturn, trace = FALSE, include.mean = FALSE)
summary(btc.garch1)
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~arma(3, 2) + garch(1, 1), data = bitcoin.logreturn, 
##     cond.dist = "std", include.mean = FALSE, trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ arma(3, 2) + garch(1, 1)
## <environment: 0x7febe63efae8>
##  [data = bitcoin.logreturn]
## 
## Conditional Distribution:
##  std 
## 
## Coefficient(s):
##       ar1        ar2        ar3        ma1        ma2      omega  
##  0.987930  -0.047911   0.052645  -1.000000   0.030285   1.039711  
##    alpha1      beta1      shape  
##  0.960721   0.745631   2.230441  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## ar1      0.98793     0.67436    1.465   0.1429    
## ar2     -0.04791     0.66602   -0.072   0.9427    
## ar3      0.05264     0.02293    2.295   0.0217 *  
## ma1     -1.00000     0.67724   -1.477   0.1398    
## ma2      0.03028     0.66627    0.045   0.9637    
## omega    1.03971     0.73248    1.419   0.1558    
## alpha1   0.96072     0.62862    1.528   0.1264    
## beta1    0.74563     0.02867   26.004   <2e-16 ***
## shape    2.23044     0.18027   12.373   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  -4809.003    normalized:  -2.63507 
## 
## Description:
##  Wed Dec 26 14:15:20 2018 by user:  
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value     
##  Jarque-Bera Test   R    Chi^2  3351.46   0           
##  Shapiro-Wilk Test  R    W      0.9080195 0           
##  Ljung-Box Test     R    Q(10)  37.95238  3.868291e-05
##  Ljung-Box Test     R    Q(15)  48.73311  1.935663e-05
##  Ljung-Box Test     R    Q(20)  51.61716  0.0001294451
##  Ljung-Box Test     R^2  Q(10)  7.820492  0.6463659   
##  Ljung-Box Test     R^2  Q(15)  15.60675  0.4086575   
##  Ljung-Box Test     R^2  Q(20)  20.57296  0.422639    
##  LM Arch Test       R    TR^2   9.994276  0.6164628   
## 
## Information Criterion Statistics:
##      AIC      BIC      SIC     HQIC 
## 5.280003 5.307173 5.279955 5.290026
summary(btc.garch2)
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~arma(3, 2) + garch(1, 1), data = bitcoin.logreturn, 
##     cond.dist = "sstd", include.mean = FALSE, trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ arma(3, 2) + garch(1, 1)
## <environment: 0x7febe4675cf0>
##  [data = bitcoin.logreturn]
## 
## Conditional Distribution:
##  sstd 
## 
## Coefficient(s):
##       ar1        ar2        ar3        ma1        ma2      omega  
##  0.986282  -0.058434   0.057549  -1.000000   0.037495   1.059835  
##    alpha1      beta1       skew      shape  
##  1.000000   0.748162   0.950560   2.219863  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## ar1      0.98628     0.59577    1.655   0.0978 .  
## ar2     -0.05843     0.58216   -0.100   0.9200    
## ar3      0.05755     0.02317    2.483   0.0130 *  
## ma1     -1.00000     0.59873   -1.670   0.0949 .  
## ma2      0.03749     0.58322    0.064   0.9487    
## omega    1.05983     0.73726    1.438   0.1506    
## alpha1   1.00000     0.64617    1.548   0.1217    
## beta1    0.74816     0.02820   26.527   <2e-16 ***
## skew     0.95056     0.01897   50.106   <2e-16 ***
## shape    2.21986     0.16825   13.194   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  -4805.75    normalized:  -2.633288 
## 
## Description:
##  Wed Dec 26 14:15:23 2018 by user:  
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value     
##  Jarque-Bera Test   R    Chi^2  3407.956  0           
##  Shapiro-Wilk Test  R    W      0.9074388 0           
##  Ljung-Box Test     R    Q(10)  38.58948  2.995275e-05
##  Ljung-Box Test     R    Q(15)  49.07985  1.700474e-05
##  Ljung-Box Test     R    Q(20)  51.90263  0.00011764  
##  Ljung-Box Test     R^2  Q(10)  7.778186  0.6504932   
##  Ljung-Box Test     R^2  Q(15)  15.5053   0.4156685   
##  Ljung-Box Test     R^2  Q(20)  20.5844   0.4219463   
##  LM Arch Test       R    TR^2   9.932961  0.621842    
## 
## Information Criterion Statistics:
##      AIC      BIC      SIC     HQIC 
## 5.277535 5.307723 5.277475 5.288670
btc.garch3 = garchFit(~garch(1,1), cond.dist = "std", data = bitcoin.logreturn, trace = FALSE, include.mean = FALSE)
btc.garch4 = garchFit(~garch(1,1), cond.dist = "sstd", data = bitcoin.logreturn, trace = FALSE, include.mean = FALSE)
summary(btc.garch3)
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~garch(1, 1), data = bitcoin.logreturn, cond.dist = "std", 
##     include.mean = FALSE, trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## <environment: 0x7febe3fc12b0>
##  [data = bitcoin.logreturn]
## 
## Conditional Distribution:
##  std 
## 
## Coefficient(s):
##   omega   alpha1    beta1    shape  
## 0.69312  0.69577  0.74137  2.37397  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## omega    0.69312     0.34260    2.023   0.0431 *  
## alpha1   0.69577     0.27471    2.533   0.0113 *  
## beta1    0.74137     0.02812   26.365   <2e-16 ***
## shape    2.37397     0.18533   12.810   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  -4832.243    normalized:  -2.647805 
## 
## Description:
##  Wed Dec 26 14:15:24 2018 by user:  
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value     
##  Jarque-Bera Test   R    Chi^2  3693.107  0           
##  Shapiro-Wilk Test  R    W      0.9077907 0           
##  Ljung-Box Test     R    Q(10)  57.00746  1.32864e-08 
##  Ljung-Box Test     R    Q(15)  76.57821  2.932977e-10
##  Ljung-Box Test     R    Q(20)  85.4443   4.58843e-10 
##  Ljung-Box Test     R^2  Q(10)  7.392301  0.6879625   
##  Ljung-Box Test     R^2  Q(15)  18.34308  0.2450587   
##  Ljung-Box Test     R^2  Q(20)  23.50332  0.2647611   
##  LM Arch Test       R    TR^2   9.92011   0.6229693   
## 
## Information Criterion Statistics:
##      AIC      BIC      SIC     HQIC 
## 5.299993 5.312068 5.299983 5.304447
summary(btc.garch4)
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~garch(1, 1), data = bitcoin.logreturn, cond.dist = "sstd", 
##     include.mean = FALSE, trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## <environment: 0x7febe6c46f20>
##  [data = bitcoin.logreturn]
## 
## Conditional Distribution:
##  sstd 
## 
## Coefficient(s):
##   omega   alpha1    beta1     skew    shape  
## 0.93802  0.98709  0.74297  0.91947  2.24705  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## omega    0.93802     0.62476    1.501    0.133    
## alpha1   0.98709     0.60047    1.644    0.100    
## beta1    0.74297     0.02793   26.599   <2e-16 ***
## skew     0.91947     0.01751   52.517   <2e-16 ***
## shape    2.24705     0.17856   12.584   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  -4822.316    normalized:  -2.642365 
## 
## Description:
##  Wed Dec 26 14:15:25 2018 by user:  
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value     
##  Jarque-Bera Test   R    Chi^2  3686.413  0           
##  Shapiro-Wilk Test  R    W      0.9076751 0           
##  Ljung-Box Test     R    Q(10)  57.1888   1.2284e-08  
##  Ljung-Box Test     R    Q(15)  76.796    2.677925e-10
##  Ljung-Box Test     R    Q(20)  85.70244  4.140708e-10
##  Ljung-Box Test     R^2  Q(10)  7.538361  0.6738286   
##  Ljung-Box Test     R^2  Q(15)  17.96327  0.2646013   
##  Ljung-Box Test     R^2  Q(20)  23.15379  0.2813062   
##  LM Arch Test       R    TR^2   10.09792  0.6073708   
## 
## Information Criterion Statistics:
##      AIC      BIC      SIC     HQIC 
## 5.290210 5.305304 5.290195 5.295778
btc.garch5 = garchFit(~arma(1,3)+garch(1,1), cond.dist = "std", data = bitcoin.logreturn, trace = FALSE, include.mean = FALSE)
btc.garch6 = garchFit(~arma(1,3)+garch(1,1), cond.dist = "sstd", data = bitcoin.logreturn, trace = FALSE, include.mean = FALSE)
summary(btc.garch5)
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~arma(1, 3) + garch(1, 1), data = bitcoin.logreturn, 
##     cond.dist = "std", include.mean = FALSE, trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ arma(1, 3) + garch(1, 1)
## <environment: 0x7febe41877e8>
##  [data = bitcoin.logreturn]
## 
## Conditional Distribution:
##  std 
## 
## Coefficient(s):
##       ar1        ma1        ma2        ma3      omega     alpha1  
##  0.992415  -1.000000  -0.024200   0.054454   1.025406   0.945721  
##     beta1      shape  
##  0.745709   2.234780  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## ar1      0.99242     0.00621  159.811   <2e-16 ***
## ma1     -1.00000     0.02341  -42.717   <2e-16 ***
## ma2     -0.02420     0.02994   -0.808   0.4189    
## ma3      0.05445     0.02160    2.521   0.0117 *  
## omega    1.02541     0.70398    1.457   0.1452    
## alpha1   0.94572     0.59961    1.577   0.1147    
## beta1    0.74571     0.02834   26.308   <2e-16 ***
## shape    2.23478     0.17725   12.608   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  -4808.799    normalized:  -2.634958 
## 
## Description:
##  Wed Dec 26 14:15:26 2018 by user:  
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value     
##  Jarque-Bera Test   R    Chi^2  3364.688  0           
##  Shapiro-Wilk Test  R    W      0.9079091 0           
##  Ljung-Box Test     R    Q(10)  36.52551  6.834724e-05
##  Ljung-Box Test     R    Q(15)  47.01264  3.66483e-05 
##  Ljung-Box Test     R    Q(20)  49.6626   0.0002474824
##  Ljung-Box Test     R^2  Q(10)  7.792423  0.6491045   
##  Ljung-Box Test     R^2  Q(15)  15.44085  0.4201527   
##  Ljung-Box Test     R^2  Q(20)  20.44682  0.4303092   
##  LM Arch Test       R    TR^2   9.93185   0.6219394   
## 
## Information Criterion Statistics:
##      AIC      BIC      SIC     HQIC 
## 5.278684 5.302834 5.278645 5.287592
summary(btc.garch6)
## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~arma(1, 3) + garch(1, 1), data = bitcoin.logreturn, 
##     cond.dist = "sstd", include.mean = FALSE, trace = FALSE) 
## 
## Mean and Variance Equation:
##  data ~ arma(1, 3) + garch(1, 1)
## <environment: 0x7febdd6dca08>
##  [data = bitcoin.logreturn]
## 
## Conditional Distribution:
##  sstd 
## 
## Coefficient(s):
##       ar1        ma1        ma2        ma3      omega     alpha1  
##  0.985875  -1.000000  -0.022856   0.059004   1.059839   1.000000  
##     beta1       skew      shape  
##  0.748338   0.950548   2.219593  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## ar1     0.985875    0.008246  119.559  < 2e-16 ***
## ma1    -1.000000    0.023904  -41.835  < 2e-16 ***
## ma2    -0.022856    0.030004   -0.762  0.44621    
## ma3     0.059004    0.021336    2.765  0.00569 ** 
## omega   1.059839    0.737041    1.438  0.15044    
## alpha1  1.000000    0.645360    1.550  0.12126    
## beta1   0.748338    0.027891   26.831  < 2e-16 ***
## skew    0.950548    0.018939   50.189  < 2e-16 ***
## shape   2.219593    0.167097   13.283  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  -4805.513    normalized:  -2.633158 
## 
## Description:
##  Wed Dec 26 14:15:28 2018 by user:  
## 
## 
## Standardised Residuals Tests:
##                                 Statistic p-Value     
##  Jarque-Bera Test   R    Chi^2  3405.352  0           
##  Shapiro-Wilk Test  R    W      0.9074002 0           
##  Ljung-Box Test     R    Q(10)  38.94565  2.595094e-05
##  Ljung-Box Test     R    Q(15)  49.25859  1.590455e-05
##  Ljung-Box Test     R    Q(20)  51.97689  0.0001147455
##  Ljung-Box Test     R^2  Q(10)  7.784134  0.6499131   
##  Ljung-Box Test     R^2  Q(15)  15.38353  0.4241596   
##  Ljung-Box Test     R^2  Q(20)  20.46078  0.4294574   
##  LM Arch Test       R    TR^2   9.923766  0.6226486   
## 
## Information Criterion Statistics:
##      AIC      BIC      SIC     HQIC 
## 5.276179 5.303348 5.276131 5.286201
#btc.garch1.res=btc.garch1@residuals    # residuals time series
#btc.garcho.res=btc.garcho@residuals    # residuals time series
#acf(btc.garch1.res); acf(btc.garch1.res**2)
#acf(btc.garcho.res); acf(btc.garcho.res**2)
par(mfrow=c(2,1))
btc.garch5.res=btc.garch5@residuals # residuals time series
acf(btc.garch5.res,lag.max = 100); acf(btc.garch5.res**2,lag.max = 100)

btc.garch5.pred=predict(btc.garch5,n.ahead=300,trace=F, plot = TRUE)

Literature Review of Univariate Garch model

Data Manipulation

NA data equals previous day

library(timeSeries)
# Equal non-existent data to previous day's
data.previous <- interpNA(data, method = "before")
# calculate log return
dataset.log.previous <- na.trim(diff(log(data.previous)))
plot.zoo(dataset.log.previous, main = "Returns with NA data Equaling Previous Day")

CCF

for (m in (2:ncol(dataset.log.previous))) {
  ccfvalues <- ccf(dataset.log.previous[,m], dataset.log.previous[,1])
  ccfvalues
  lag2.plot(dataset.log.previous[,1], dataset.log.previous[,m],10)
}

Remove row with NA

library(timeSeries)
# remove rows with non-existent data
data.complete <- removeNA(data)
# calculate log return
dataset.log.complete <- removeNA(diff(log(data.complete)))
plot.zoo(dataset.log.complete, main = "Returns with NA Data Removed")

fit model

\[X_t = \beta'z_t+y_t^* = \beta_1 Fed_{t-1} + \beta_2 EURUSD_{t-1} + \beta_3 DJI_{t-1} + \beta_4 Gold_{t-1} + \beta_5 Oil_{t-1} + \beta_6 RUT_{t-1} + \beta_7 GSPC_{t-1} + \beta_8 EEM_{t-1}\]

\[\phi(B) y_t^* = \theta(B) y_t\] or

\[y_t^* = \beta_0 + \beta_1 y_{t-1} + \beta_2 y_{t-2}\]

\[y_t = \sigma_t \epsilon_t\]

\[\sigma_t^2 = \alpha_0 + \sum_{j=1}^m{\alpha_j}y_{t-j}^2 + \sum_{j=1}^s{\beta_j}\sigma_{t-j}^2=\alpha_0+\alpha y_{t-1}^2 + \beta \sigma_{t-1}^2\]

dlp <- as.data.frame(dataset.log.previous) 

library(dplyr); library(knitr)

modeldata <- mutate(dlp, z1 = lag(Federal.Rate,7),z2 = lag(EURUSD,2), z3 = lag(DJI,6),z4 = lag(Gold,5),z5 = lag(Crude.Oil,7),z6 = lag(RUT,6),z7 = lag(GSPC,6),z8 = lag(EEM,9))
modeldata1 <- modeldata %>% select(Adj.Close,z1, z2, z3, z4, z5, z6, z7, z8)

modelfit <- lm(Adj.Close ~ ., data = modeldata1)
summary(modelfit)
res <- modelfit$residuals
par(mfrow=c(1,2))
qqnorm(res);qqline(res);acf(res)
auto.arima(res, seasonal = FALSE, ic = "bic")
resarma <- sarima(res, 3, 0, 2)
resss <- resid(resarma$fit)
par(mfrow=c(2,1))
acf(resss);acf(resss**2,lag.max = 100)
source("garchAuto.R")
fit1 = garchAuto(res, cores=8, trace=TRUE) ## can only work in MAC
fit1
res.garch5 = garchFit(~arma(1,3)+garch(1,1), cond.dist = "std", data = res, trace = FALSE, include.mean = FALSE)

par(mfrow=c(2,1))
res.garch5.res=res.garch5@residuals # residuals time series
acf(res.garch5.res,lag.max = 100); acf(res.garch5.res**2,lag.max = 100)